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Abstract 

(N. 

We describe a new mechanism for pulsations in evolved stars: relaxation oscillations driven by a coupling 
between the luminosity-dependent mass-loss rate and the H fuel abundance in a nuclear-burning shell. When 
mass loss is 

included, the outward flow of matter can modulate the flow of fuel into the shell when the stellar 
luminosity is close to the Eddington luminosity LEdd- When the luminosity drops below Z/Edd, the mass 
outflow declines and the shell is re-supplied with fuel. This process can be repetitive. We demonstrate 
the existence of such oscillations and discuss the dependence of the results on the stellar parameters. In 
particular, we show that the oscillation period scales specifically with the mass of the H-burning relaxation 
shell (HBRS), defined as the part of the H-burning shell above the minimum radius at which the luminosity 
from below first exceeds the Eddington threshold at the onset of the mass loss phase. For a stellar mass 
£N| M* ~ 0.7M Q , luminosity L„ ~ 10 4 L Q , and mass loss rate \M\ ~ 1O _5 M0 yr _1 , the oscillations have a 

recurrence time ~ 1400 years ~ 57Tf sm , where Tf sm is the timescale for modulation of the fuel supply in the 
HBRS by the varying mass-loss rate. This period agrees very well with the ~ 1400-year period inferred 
for the spacings between the shells surrounding some planetary nebulae. We also find the half-width of the 
luminosity peak to be ~ 0.39 times the oscillation period; for comparison, the observational shell thickness 
of ~ 1000 AU corresponds to ~ 0.36 of the spacing between pulses. We find oscillations only for models 
in which the luminosity of the relaxation shell is ~ 10-15% of the total stellar luminosity and for which 
energy generation occurs through the pp chain. We suggest this mechanism as a natural explanation for the 
circumnebular shells surrounding some planetary nebulae, which appear only at the end of the AGB phase. 

^ ' 1. Introduction 

Asymptotic giant branch (AGB) stars exhibit pulsations on a variety of timescales, ranging from long 
periods associated with thermal pulses of the He-burning shells (r ~ 10 5 years; c/., Iben and Renzini 1983, 
Schonberner 1983, Mazzitelli and D'Antona 1986, Iben 1991, Dorman el ah, 1993, Vassiliadis and Wood 
1993) to far shorter oscillations associated with acoustic modes of the star (r ~ 1 year). Recently, a 
new class of pulsational behavior has been inferred for such stars. The evidence takes the form of nested 
circumnebular shells around a number of planetary nebulae (c/., Mauron and Huggins 2000; Terzian and 
Hajian 2000; Balick, Wilson, and Hajian 2001). These appear to have been ejected in pulses during the AGB 
phase or succeeding post- AGB but pre-planetary-nebula (pre-PN) phases of stellar evolution. For example, 
surrounding the Cat's Eye Nebula (NGC 6543), Balick et al. (2001) find rings that appear to have been 
ejected quasi-periodically, with the estimated time between successive pulses ranging from ~ 1000 to ~ 1900 
years; the mean spacing for eight successive shells is approximately 1400 years. The authors take the rings 
to be ~ 1000 AU thick, indicating that the ejection mechanism operates during 0.36 of the time between 
successive pulses. Each shell is estimated to contain ~ 0.01M Q of material, and a total mass ~ 0.1M© 
appears to have been ejected in the process that produced the shells. 

To date, four models have been proposed to explain these enigmatic shells. Soker and Harpez (1988) sug- 
gested a binary model, in which the secondary follows a highly elliptical orbit. As the secondary approaches 
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pcriastron, it interrupts the mass loss from the primary, leading to the formation of a shell. The intershcll 
timescale is then set by the orbital period r or b of the binary. Mastrodemos and Morris (1998, 1999) proposed 
a different binary model, in which the shells are the result of spiral shocks in the AGB wind. In these models 
the intershell timescale is also r or b. Simis, Icke, and Dominik (2001) proposed that the shells are produced 
by quasi-periodic variations in dust formation rates in dust-driven AGB winds. They argue that non- linear 
coupling between pulsation of the star, dust formation in the atmosphere, and radiation-driving in the wind 
produces shells with a quasi-regular period ~ 10 3 years. Finally, Garcia-Segura et al. (2002) have proposed 
an MHD model in which a dynamo operating in the star sets the period of the oscillations, and magnetic 
pressure in the wind creates the shells. It is difficult at present to evaluate the validity of these models. The 
eccentric binary model seems unlikely to account for the several PN discovered with multiple shells; it is not 
clear why a ~ 10 3 year period should be favored; and it is hard to understand how the observed variation 
from ~ 1000 to 2000 years in shell spacings could be produced by a perfectly periodic binary model. The 
MHD model may be based on sound physics, but existing AGB dynamo models (Blackman et al. 2001) 
predict far shorter dynamo periods. 

A number of indications seem to us to point instead to fuel-supply-limited relaxation oscillations as the 
explanation for the observed circumnebular mass shells: 

(1) At high luminosities, stars lose mass at exceptionally high rates. Indeed, the brightest precursors 
of the planetary nebulae may have L» ^ L^dd, driving the so-called "superwind" that terminates evolution 
up the AGB. However, the effect of the stellar wind upon the internal structure of a star has not, to our 
knowledge, been included self-consistently in previous calculations. Models of the wind assume that the 
stellar interior is dynamically unaffected by the mass flux through the surface layers. Conversely, models of 
stellar evolution that include wind mass loss simply strip mass away from the stellar surface. Because the 
rate of mass loss increases dramatically as a star climbs the AGB (c/., Vassiliadis and Wood 1993), any effect 
due to the wind will be accentuated near the end of the AGB phase. As the observed circumnebular shells 
were evidently produced during the final stages of evolution just prior to the ejection of a planetary nebula, 
some effect of mass loss upon the internal stellar structure thus seems a natural candidate to explain the 
pulsed ejection of the circumnebular shells. 

(2) Comparison of the rate at which H burning consumes mass with the rate of mass loss from the stellar 
surface also is suggestive. Let Q' be the amount of energy released per gram of H consumed. This quantity 
is (c/., Clayton 1968, p. 287 ff) 

Q> = (4mH ' mHc)c2 = 6.398 x 10 18 erg g ~\ (1) 
4toh 

If X is the mass fraction of H in the unburned material, M s the rate at which an H-burning shell consumes 
mass, and L s the resulting luminosity produced by this shell, then 

M.X = -^. (2) 

With L s ~ 10 4 Lq and X ~ 0.7, the rate at which mass must be consumed by a H-burning shell in order to 
provide the stellar luminosity is thus \M S \ <~ 10~ 7 M Q yr _1 . For comparison, the observed rates of mass loss 
range from 5 x 10~ n (\M\/M & yr _1 ) ^ 2 x 10~ 6 for the central stars of the planetary nebulae (Pcrinotto 
1989) to - 10~ 4 M Q yr" 1 for long-period variables (c/., Cox et al. 2000, p. 415). We note that the OH/IR 
stars which have the highest mass loss rates may have L* ^ i^Edd- Thus, at most stages in an AGB star's 
evolution, mass loss from the surface is comparable to or higher than mass consumption at the H-burning 
shell. Consequently, there is a competition between the need to supply fuel downward into the H-burning 
shell and upward into the stellar wind. 

We are thus led to suggest that an interaction between radiation-driven mass loss from the stellar surface 
and mass consumption by H-burning in a highly evolved AGB star may produce a relaxation oscillation that 
ejects shells of matter quasi-periodically. Near the tip of the AGB, a star has a luminosity ~ 10 4 L Q , and 
it contains a dense, degenerate C/O core with mass ~ 0.6M Q and radius ~ 0.02i? Q ~ 10 9 cm, surrounded by 
a tenuous H/He envelope extending outward to radius i?» ~ 4 x 10 13 cm ~ 500i? Q (c/., Ibcn 1987). Suppose 
that such a star is initially in a steady state. We define the H-burning relaxation shell (HBRS) as a thin 
outer part of the H-burning region above the smallest radius for which the luminosity from below can exceed 
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the Eddington luminosity. This region can be quite a small fraction of the total H-burning shell, but is the 
region which is important for relaxation oscillations that we now discuss. 

We assume that the HBRS consumes mass at a rate \M S \. Suppose also that the star undergoes radiation- 
driven mass loss from the stellar surface at a steady rate |M|. A perturbation that increases the temperature 
of the HBRS, leading to a thermonuclear runaway, greatly increases the luminosity. Because the surface 
mass-loss rate is proportional to the luminosity (c/., §2), the increased luminosity increases the rate of mass 
loss from the surface. Conservation of mass, however, requires that the mass flowing out through the surface 
must be replaced from below; that is, mass loss is not simply a stripping away of the outer layers of the 
star. Some form of gradual readjustment must occur along the AGB and into the post-AGB and pre-PN 
phases. The increased mass-loss rate thus at least impedes the supply of fuel to the HBRS, and it may even 
temporarily interrupt that supply. In any event, the increased surface mass loss has the effect of reducing 
the luminosity produced by the HBRS, which in turn reduces the rate of radiation-driven mass loss from 
the stellar surface, allowing the cycle to repeat itself. This is a previously unexplored coupling between the 
HBRS luminosity, the mass-loss rate, and the fuel supply available in the HBRS. Such fuel-supply-limited 
relaxation oscillations are familiar in other contexts, such as the pulsation of a fuel-starved flame from a 
guttering candle. 

In the present paper, we consider in more detail the hypothesis that the source of the circumncbular 
shells is just such a fuel-supply-limited relaxation oscillation in a star near the end of its AGB phase. In 
§2, we first review briefly the relation between the rate of mass loss M and the stellar luminosity L* for 
the optically thick, radiation-driven winds appropriate to high-luminosity stars. In §3, we next collect 
together the fundamental Eulerian equations that describe the hydrodynamic interaction between matter 
and radiation in a spherically symmetric star with a spherically symmetric mass outflow. In §4 we develop 
a simplified model that accounts qualitatively for the main features of the relaxation oscillations described 
above. We emphasize that this is a simplified model, intended to illustrate the essential physics of these 
oscillations. A more complete calculation, beyond the scope of this paper, would require a detailed stellar 
evolution model including self-consistent mass loss. We describe the results of numerical calculations using 
our simplified model in §5, and we conclude in §6 with a summary. 

2. Mass Loss on the AGB 

From his observations of mass loss from luminous K and M giants and supergiants, Reimers (1975a, b) 
proposed an empirical formula to express the dependence of the rate of mass loss upon the stellar parameters: 

where r\ is a parameter of order unity. This expression has been used extensively, both in analyzing ob- 
servations and in theoretical studies. For example, for an AGB star with M* ~ Mq, L* ~ IO^Lq, and 
i?* <~ 5OOi?0, this expression yields the result \M\ ~ 2 x 10~ 6 M Q yr" 1 . However, as originally pointed 
out by Renzini (1981; see also Iben and Renzini 1983), this equation appears to underestimate the rate of 
mass loss in the very high luminosity, thermally pulsing AGB phase, requiring a postulated "superwind" to 
provide the necessary mass-loss rates. 

Bowcn (1988) and Bowen and Willson (1991) have described hydrodynamic calculations of the atmo- 
spheres of long-period, Mira variables - pulsating AGB stars of low to intermediate masses - that lead to 
significant rates of mass loss. In models that include dust formation, radiation pressure acting on the dust 
grains produces a mass loss extending up to M ~ 10 _6 M Q yr -1 for a star with = 1.2M Q , — 270i? Q , 
and = 5315-Lq. Bowen suggests this as the "superwind" that appears to terminate AGB evolution. 
Willson (2000) has recently re-interpreted the empirical Reimers mass-loss rate (3) in terms of these results. 

In their computations of the evolution of intermediate-mass stars, Vassiliadis and Wood (1993) included 
an empirical mass-loss rate that reproduces the very large values of M observed during the AGB phase of 
evolution. It reaches an essentially constant "superwind" rate that is within a factor of two of the value 

|M| = £ ~ 1O-M yr- x ^ (4) 
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appropriate for a radiation-driven wind, where v w ~ 15 km s _1 is the asymptotic wind speed far from the 
central star. For stars with M ^ 2.5M , the "superwind" mass-loss rate is achieved only during the last few 
thermal pulses, and the vast majority of the mass loss occurs in these episodes. 

Koesterke, Dreizler, and Rauch (1998) and Koesterke and Werner (1998) measured the mass-loss rates 
for four PG 1159 stars - post-AGB stars that are among the immediate precursors of the white dwarfs - 
finding —8 & log(|M|/M© yr _1 ) ^ —7, consistent with the theory of radiation-driven winds. 

3. Fundamental Equations 

The fundamental Eulerian equations that describe the structure and time variation of a spherically 
symmetric star with spherically symmetric mass loss are (c/., Landau and Lifshitz 1959; Hansen and Kawaler 
1994) the continuity equation, 



dp 1 d(r 2 pv) 
dt r 2 dr 

where p is the density and v is the radial velocity; the definition of the mass M r interior to radius r, 

dM, 



dr 

the equation of conservation of radial momentum, 



= Anr A p: (56) 



dP fdv dv\ GM r 

+ p \^i + v ^r =-P— 5-' ( 5c ) 



dr \dt dr J 

where P = P g + P r is the sum of the gas pressure P g = pksT/ pH and radiation pressure P r — |a r T 4 , and 
where T is the temperature; and the heat-flow equation 



dT 

dr 



3 Kp L r . . . 

^j- — m radiative equilibrium; 



4ac T* 4irr z 

w TdP . ru . (M) 

Vad^"^! m convectivc equilibrium. 
P dr 



Here k is the opacity of the stellar matter, and L r is the luminosity flowing out through a sphere of radius 
r. Note that the first form of equation (5d) assumes that energy transport is dominated by radiation flow 
rather than convection, while the second form is that appropriate for convective energy transport. We also 
require the equation of energy conservation, 

9L r A 2 ( _ dS\ ._ , 

— = 4^r 2 p e - T— - Tv— , (5e) 



dr \ dt dr 

where e is the thermonuclear energy generation rate, and S is the entropy. In addition, we need the equation 
for the time rate of change of the H mass fraction X due to nuclear burning in the HBRS: 



dX dX _ e 
~dt +V ~dr~ ~ ~Q 



v— = -7v. (5/) 



where Q' is given by equation (1). Note that if we multiply equation (5f) by dM r and integrate over the 
HBRS, we recover equation (2), with M s = 4irr 2 p s v s . 

4. A Simplified Model 

Integrating equation (5b) from the center of the star out to some fixed radius r and computing the time 
rate of change of this quantity gives 

^ r _dM L= f\ nr 2^P dr = _ Al:r 2 pv = - F ^ t ) i ( 6 ) 
dt In dt 
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where we use equation (5a) to eliminate dp/dt, and the last equivalence defines the local mass-flow rate F. At 
the stellar surface, M r — > M = —F(R tl [t],t), a quantity that depends only upon time. If the star experiences 
steady-state mass loss, then equation (5a) shows that the mass flux F(r,t) is a constant, independent of 
position as well as time. In general, F may depend upon both position and time, but if the changes in the 
star occur quasi-statically (as is the case in stellar evolution calculations), then we can expect the mass flux 
through the star to remain close to its (position-independent) value throughout the entire stellar envelope. 

From a scale analysis of the momentum equation (5c), we find that the ratio of the velocity-dependent 
terms to the remaining terms is of order (v/c s ) 2 , where c s is the sound speed. For a mass-loss rate of 
10 _5 M Q yr _1 , and assuming the mass flow rate F = Anr 2 pv to be independent of position in the stellar 
envelope, the flow velocity at the HBRS is v ~ 0.5 cm s" 1 . In comparison, the sound speed at the shell is 
c s ~ 4 x 10 7 cm _1 . Thus, at the HBRS, the velocity-dependent terms in equation (5c) are only ~ 10~ 16 
times the other terms in this equation. These terms also appear to be negligible near the photosphere, and 
for this reason, we neglect them in the momentum equation. This approximation is equivalent to assuming 
that the variations of interest to us here all take place on timescales long in comparison with the hydrostatic 
readjustment timescale. For this reason, the hydrodynamic timescale Th y d ~ R*/v ~ 500i?©/15 km s _1 ~ 1 
year plays no further role in our development of the equations that describe this model. 

Dropping these w-dependent terms reduces the momentum equation to the usual equation of hydrostatic 
equilibrium. Separating out the radiation pressure P r , we find that the equation of hydrostatic equilibrium 
for the gas pressure P g can be written as 



dP g _ pGM r 
dr r 2 



4ncGM r 



(7) 



In the envelope of the star we can reasonably assume L r w L* = constant and M r w M* = constant. With 
these approximations, equations (7) and (5d) can be written in the forms 



dP g pGM* 
dr r 2 



AncGM* 



= -o— (1-A), (8a) 



where the final equivalence defines the parameter A, and we have assumed radiative energy transport, at 
least near the HBRS: 

dT = 3 up L* 

dr 4a r cT 3 47rr 2 ' 1 ' 

The quantity A in equation (8a) has a simple physical interpretation - it is a dimensionless form of the stellar 
luminosity: 

X = kL * = L * (q) 
A-kcGM* L Ed d' 

where L Edd = AttcGM^/k is the Eddington luminosity. 

The assumption that radiative transport dominates in the stellar envelope is not clearly justified a 
priori, as the highly distended surface layers of AGB stars are well known to contain deep convection zones. 
We nevertheless believe that equation (8b) is appropriate here for at least two reasons: (1) The surface 
convection zone in an AGB star does not extend all the way down to the HBRS. We are most concerned 
with conditions near this shell, and it is well known that the temperature profile below the convection zone 
in the surface layers of a star rapidly approaches the "radiative zero" approximation we use below. (2) It is 
not clear whether the relaxation oscillations we wish to study actually occur on the AGB or whether they 
occur at some later high-luminosity phase, where the surface convection zone may be shallower and equation 
(8b) correspondingly better justified. 

Dividing equation (8a) by equation (8b), assuming a Kramers' law opacity of the form* 

K = K*P%P-\ (10) 



* The more familiar form of this equation is k = K pT 3 5 . The form given in equation (10) is more 
convenient for our present purposes, however, and the two forms are equivalent if a = 1 and b = 1.125. 
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where k„, a, and 6 are constants, and employing the conventional definition 



p g /p r = a/{i-(i) (ii) 

enables us to write the equation of hydrostatic equilibrium in a form that can be integrated directly, if wc 
assume (3 to be independent of position in the star. Using the "radiative zero" surface boundary condition, 
P g = at r = i?*, we obtain 

(6 - a + 1) 

Note that equation (12) includes the full effect of the Kramers' law opacity, subject only to the approximation 
that (3 is independent of position. 

Using equation (10) for k in equation (8b), expressing p in terms of P g and T with the ideal gas law, 
and again using equation (11) with (3 assumed independent of position, we can also integrate equation (8b), 
obtaining for the temperature distribution in the stellar envelope 

{b-a+ |) GM*nH (I 1 



T (b-a + 1) k [r i?J' (13) 

where we assume T = at r = R*, consistent with the approximation made in deriving equation (12). 
Equation (13) expresses the temperature at position r in terms of constants and stellar properties, including 
the stellar radius 

Equations (12) and (13) express the conditions of hydrostatic and thermal equilibrium in the stellar 
envelope. Using these equations, we can also obtain an expression for the mass contained in the stellar 
envelope above the HBRS at radius r s : 



AM., 



r 



A7rpr 2 dr 



k (1-/3) 3 * r 



b — a 



b-a+1 



GM^H 



-, 3 



kR* 



IvJr, (« 



s ds, 



(14) 



where s = r/R*. The integral in equation (14) can be evaluated analytically and yields 



/ ( 



i-1 

S 



s ds = — In — 3 

-ft* 



1- ^ 



i - 4e- 



(15) 



Note that the factors of i?* cancel in the constants in equation (14), so that the envelope mass depends 
only upon (3 (and thus, from equation [12] upon A), M*, and the fractional radius r s /R*. Consequently, as 
the stellar luminosity increases (A increases and (3 accordingly decreases), the parameter r s /R* must vary 
in order to keep AM env fixed. From the ratio of the quantity AM env in the perturbed state to the quantity 
AMcnl = AM onv in the unperturbed state, together with equation (12), we obtain 



J_ -ln 5 ~3(l-.s) + f(l-s 2 )-l(l- S 3 ) 
A(°) " - Inao - 3(1 - so) + §(!-*§)- |(1 



(16) 



where now s = r s /R*, and so is the value of s in the unperturbed state. If r s ~ 10 10 cm and _R* ~ 10 12 
cm - e.g., corresponding to a star leaving the AGB - then s ~ 10~ 2 . If we denote the temperature of the 



HBRS by T s , we can use equation (13) to express it in the form of a ratio to the temperature T 6 
unperturbed state: 



(o) • 



y = 



i-i 

s 



,(o) 



1-1 

so 



in the 



(17) 
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Note that equations (16) and (17) provide a parametric relation for the quantity X/X^ in terms of the 
dimensionlcss temperature y and the constant sq. 

Next, consider the energy equation, using the form of TdS that includes both radiation and gas pressure. 
The general expression can be written in the form (c/., Cox and Giuli 1968, p. 206 ff) 



T 



<it -r, \)d P 



TdS = C v 

Here the quantity (r 3 — 1) is the thermodynamic derivative 

(T 3 - 1) 



(18) 



2 (4 - 3/3) 



<91nT\ 

dlnp J s ~ 3 (8 - 7/3) ' 

and the last equality gives the explicit expression for a mixture of gas plus radiation. Using equation (18), 
the energy equation (5e) becomes 



dr 



= 4wr p 



{- 



dT T dp 
0t ~7 ( 3 ^ 



-vC v 



dT T dp 

"a V 1 3 - tj"^- 

or p or 



We now wish to integrate this equation over the stellar envelope. If we assume 

e = e X P T v , 



(19) 



(20) 



we can integrate the first term on the right-hand side of equation (19) to obtain the luminosity generated 
by the HBRS, 



edM r 




(21a) 



Here the subscript s denotes conditions at the HBRS at some time t, and the superscript (0) denotes the 
conditions in the unperturbed state. The temperature exponent v in the H-burning reaction rate has the 
value v ~ 2 — 6 for the pp chain, which is appropriate for the relatively low shell-burning temperatures we 
consider here, while it has the value v ~ 10 — 20 for the CNO cycle (c/., Hansen and Kawaler 1994, p. 238, 
241). From equations (11) and (12) we can also obtain the relative density distribution in the envelope, 



Ps 



AW 
A 



T, 



T, 



(0) 



(22) 



and this can be used to replace p s /p^ in equation (21a). 

In the term involving time derivatives in equation (19), we use the continuity equation to eliminate 
dp/ dt. The resulting integrand involves the spatial derivative of the mass-flow rate F. We assume that F 
is independent of position as argued previously and accordingly neglect this term. The remaining integral 
involving a time derivative can be written in the approximate form 



R, flrp ,rp 

C v —dM r w C^AM cnv -^, 



dt 



dt 



(216) 



where Cy is some average heat capacity for the envelope, AM env is the envelope mass defined in equation 
(14), and T s is again the temperature of the HBRS. The remaining terms in equation (19) represent heat 
advection. We estimate them to be only ~2x 10~ 3 as large as other terms in this equation, and we have 
neglected them. 

Collecting equations (21) together, we can thus write the integral of equation (19) in the form 



v+3 



T N 



dt Vt s (0) 



(23) 



where L\ denotes the luminosity originating interior to the HBRS, and 
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(24) 



is the thermal relaxation timescale. In evaluating r t h we have taken AM onv ~ 0.1M Q ; T ~ 10 6 K, which 
may be typical of a fuel-starved HBRS; and L* ~ 10 4 L Q . 

In the final equation (5f), we first carry out a Galilean transformation to a frame moving with the HBRS. 
For any quantity f(r, t), if we define a new position coordinate to be r' = r — v s t, where v s is the propagation 
speed of the HBRS, the general relations among the partial derivatives in the old and new frames are 



df_ 

dr' 



df 
dr 



and 



dt 



dt 



' dr 



If we apply the relations (25) we can express equation (5f) in the frame moving with the HBRS: 



dX 
~dt 



+ (v - v s ) 



dX 
dr' 



e(r',t) 

Q' : 



where Q' is given by equation (1). If we define the position of the HBRS by the condition 

X(r s ,t) = ^X , 

where X ~ 0.7 is the primordial H abundance, then at the location of this shell, we have dX/dt\ Ts 
Evaluating equation (26) at this location thus gives the shell speed: 



[v{r s )-v s \ — 



Q' 



(25) 

(26) 

(27) 
= 0. 

(28) 



where v(r s ) is the speed of the flow, as measured in the rest frame of the star, at the position r s of the 
HBRS, and v s is the desired propagation speed of the shell. If we use equation (28) to describe the presumed 
unperturbed steady-flow condition and subtract this from equation (26), we obtain for the perturbed flow 
the result 



dX 
~dt 



+ (v- v s ) 



dX 
dr 7 



(#-«<°») 



dX<® 



dr' 



[e(r',f)-e«V,t)] 

Q' 



(29) 



We now need to work out the integrals of the various terms in this equation. The first term can be 
written in the form 



M, 



dX 

~dT 



dM r , = AM n 



dX 
~dt 



AM„ 



dX s 
dt 



(30a) 



This equation defines the spatial average of the time derivative, which we approximate as the time derivative 
of the H abundance in the HBRS. Note that X is constant in time and space except in the very thin 
HBRS. Accordingly the mass AM nuc in HBRS regions is relatively small; we use the approximate value 
AM nuc ~5x 1O~ 4 M for this mass. (This value is actually consistent with the final H-burning shell mass 
found by Iben (1984), Iben and MacDonald (1986), and Herwig et al. (1999) in late AGB or immediate 
post-AGB phases of stellar evolution. However since our model also requires a significant luminosity internal 
to the HBRS, L\, we do not require that the HBRS is the entire H-burning shell. Some arbitrary ratio of H 
and He could be burning below the HBRS to provide L\.) 

In the terms in equation (29) that involve spatial gradients of X, we write dM r i = 4w(r') 2 p(r' , t)dr' and 
use equation (6) to replace v with the mass-flow rate F, which we assume to be independent of position. 
The resulting integral has the form 



M, 



(v 



,dX 
Vs) d? 



(v 



(0) 



,(o) 



dr' 



dM r , = [(F - F(°>) - (F s - fM)]X , 



(306) 
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where F s = 4irr 2 pv s , and Xq is again the primordial H abundance. The integral of the terms on the 
right-hand side of equation (29) just gives the difference in the HBRS luminosities in the perturbed and 
unperturbed states. Collecting these results together, and for simplicity neglecting F s — , we write the 
integral of the X equation in the approximate form 

dX 



AM„ 



dt 



+ (F-F^)X = - 



Q' 



(31) 



The system of equations describing the presumed relaxation oscillations now consists of (23) and (31), 
with the dimensionless luminosity A expressed in terms of the dimensionless shell temperature y by the 
parametric equations (16) and (17). For given values of y and so, we solve equation (17) to obtain s, and 
we use these values in equation (16) to obtain \/\(°\ For convenience, we collect the resulting equations 
together here, defining x = Xg/X^: 



L, = Li + L ( , 0) 



/AW 



xy 



,(0) dy 
L * T th-rr, 
dt 



(0) 



{ F- F M)X = -^ r 



ft) 



xy 



v+3 



- 1 



- 1 



1 

SO 



(32a) 

(326) 
(32c) 



and 



A 

AW 



lns-3(l-s) + f(l-s 2 )-±(l-s 3 ) 



- In s - 3(1 - s ) + |(1 - s§) - |(1 - s 3 ) 
For radiation-driven mass loss, we assume the simple form given by equation (4), which we write as 



(32d) 



M = -ah*. 



(33) 



Combining (6), (9), and (33), we obtain 



M 



F 



M(0) F(°) L 



(0) 



A 



(34) 



Equation (34) eliminates equation (6) and expresses the mass flux F in terms of A. We can now put the 
equations in fully dimensionless form. Dividing equation (32a) by we obtain 



A 



= A + (1-A) 



A 



xy 



y+3 



7th 



dy 
dt' 



(35a) 



where A = Li/L^. Note that, in a truly steady state with dy/dt — 0, we have x = 1, y = 1, and 



A = A(°) , obtaining from equation (23) Lf ^ e 
Xs°^ AM nuc and use equation (34) to obtain 

dx 



dt 



+ 



(0) 
1 

Tf'sm 



■L\. Similarly, we divide equation (32b) by the quantity 



1 

7"cx 



ft) 



xy 



v+3 



- 1 



In equations (35), r t h is given by equation (24). We define 



Tf'sm — 



X s (0) AM ni 



X s (0) AM nuc 



25 years, 



(356) 



(36a) 



X Q \M\ 

as the timescale on which the fuel supply is modulated by mass loss; i.e., the time required for mass loss 
to strip material out of the HBRS. In evaluating Tf sm we assume the mass in the HBRS to be AM nuc ~ 
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5 x 10 4 Mq, as noted above, and we assume M <~ — 10 5 M Q yr 1 to be the mass-loss rate. We define the 
fuel-exhaustion timescale similarly: 

Tcx = Q'X°^Mnuc ^ x 1Q 4 yearg (366) 

An important constraint is provided by the fact that the H abundance can never exceed the primordial 
value A nor become negative. From equation (27), this correspondingly constrains our dimensionless 
parameter x = X s /X^ to the range 

< x < {Xo/^Xo) = 2. (37) 

Equations (35a) and (35b), together with the constraint equation (37), the parametric relations (32c) and 
(32d), and the definitions of the characteristic timescales (24), (36a), and (36b), constitute the system 
of equations that characterize this relaxation-oscillator model. In these equations, the quantity A*- -* is a 
dimensionless number presumably just slightly less than unity (i.e., the luminosity is just slightly less than 
the Eddington luminosity). 

The mathematical nature of the governing system of nonlinear differential equations (35a, b) is better 
seen by writing them in the fully nondimensional form 

^=A+(l-A)xh(y)-g(y), (38a) 
^ = -7ib(2/)-l]-72[*%)-l], (386) 



where 



/x A yi v +3) T th 7" th , Q( .s 

■ 9(y) "AW %) "7(^P 71 72 = (39) 



cx 



and where t is now a nondimensional time, scaled by the thermal time scale r t h- Note that g(y) = A/A' ) 
is a nonlinear function of y defined by equations (32c) and (32d). Note also that g(l) = h(\) = 1 and 
that equations (38) have the equilibrium solution x = y = 1. A linear stability analysis, presented in the 
Appendix, shows that this equilibrium is unstable to growing oscillations for a range of values of the stellar 
parameters. 

5. Numerical Solutions of the Model Equations 

We have carried out a number of numerical solutions of the nonlinear system (38) in order to explore the 
behavior of the relaxation oscillations and to determine their dependence on the various parameters in the 
model. As a specific example, consider first a benchmark model with M» = 0.7M Q , L„ = 10 4 £©, i?» = 10 13 
cm, AM onv = O.1M , AM nuc = 5x 10~ 4 M Q , and T s (0) = 10 6 K. For these parameters, equation (13) and 
the definition so = ri^/i?* yield so = 0.03. We adopt the values v = 3 and A = L\/L^ = 0.9 in this model. 
These values yield Tth = 20 years, Tf sm = 25 years, and t cx = 1.75 x 10 4 years. For these parameter choices, 
the numerical solution of equations (38) yields sustained relaxation oscillations with a period P w 72r t h ~ 
1400 years, as shown in Fig. 1. During the initial rise in temperature and luminosity, the H abundance X 
drops rapidly, reaching X = in ~ 39T t h ~ 780 years. The resulting relaxation oscillations reach a peak 
temperature T s w 6 x 10 7 K and a peak luminosity L* s=s 3.28l1°^. The half-width of a pulse, defined as the 
time during which the luminosity exceeds half its peak value, is At ~ 28T t h ~ 560 years = 0.39P. We define 
the "duty cycle" W of a pulse as the time required for the shell temperature to relax from its peak value 
to Ts°^ (i.e., for y to return to the value y = 1). This is approximately the time required for the star to 
radiate away the heat generated in the thermonuclear shell flash. A simple analytic argument based on our 
dimensionless model equations gives the duty cycle as W ~ 26T t h, while our explicit numerical calculation 
yields W = 39.4r th = 0.55P. 
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The H abundance in the shell remains vanishingly small until the luminosity - and the accompanying 
mass-loss rate - has declined sufficiently for the high-temperature shell again to overtake the H fuel supply. 
The recovery between successive pulses is governed by the time required for the H abundance - and the 
associated energy-generation rate - to recover to, and then exceed, the equilibrium rate. The time between 
pulses is P ~ 72r t h — 1400 years, closely comparable to the timescale between the observed circumnebular 
shells in NGC 6543. We obtain similar sustained relaxation oscillations starting from shell temperatures 
T s cither slightly less than or slightly greater than t]°\ where Tj ^ is the presumed steady-state burning 
temperature; in this sense, these results are robust. 

We have investigated the numerical accuracy of these results in two ways. Calculations with timesteps 
selected to be a factor of two larger and smaller than the values used in the calculations we report here give 
values of P within ±0.2%, values of W within ±1.8%, and values of the peak temperature within ±2.0%. 
Calculations for this same model obtained using two different algorithms yield values of P that agree to 
within 20%, values of W within 20%, and values of the peak temperature within about 3%. 

We have explored the variations in the properties of these relaxation oscillations for the range of model 
parameters listed in Table 1 below. Model F is the illustrative benchmark case discussed above. The 
numerical results are summarized in Table 2 below. Depending upon the specific choice of parameters, we 
obtain oscillations with a range of periods, duty cycles, and peak temperatures and luminosities. Note that 
the peak temperatures for all except model F exceed 10 8 K, violating our simplifying model assumption that 
energy production in the shell is due solely to the H-burning pp reactions. 



Table 1. Model Parameters* 



Model 


—M 


T (0) 


so 


V 


A 


Tth 


Tfsm 


Tex. 


A 


10~ 6 


10 7 


0.003 


3 


0.9 


200 


500 


3.50 x 10 4 


B 


10" 5 


10 7 


0.003 


3 


0.9 


200 


50 


3.50 x 10 4 


C 


io- B 


10 y 


0.03 


3 


0.9 


20 


500 


3.50 x 10 4 


D 


10" 5 


10 6 


0.03 


3 


0.9 


20 


50 


3.50 x 10 4 


E 


10~ 4 


10 a 


0.03 


3 


0.9 


2 


5 


3.50 x 10 4 


F 


io- & 


10* 


0.03 


3 


0.9 


20 


25 


1.75 x 10 4 



* The mass- loss rate M is in M yr 1 , is in K, and all timescales are in years. 



Table 2. Model Results* 



Model 


P 




L max/£(0) 


W 


A 


2.28 x 10 4 


31.7 


1.86 


0.45 


B 


1.38 x 10 4 


30.3 


1.85 


0.67 


C 


9.12 x 10 3 


171.7 


3.85 


0.16 


D 


2.47 x 10 3 


140.9 


3.73 


0.55 


E 


9.29 x 10 2 


1,409. 


5.05 


0.92 


F 


1.44 x 10 3 


62.2 


3.28 


0.55 



* The oscillation period P is in years, and the "duty cycle" W of the pulse is defined as the fraction of 
the pulse period required for the shell temperature T s to decline from its maximum value T s max to the 
mean steady-state temperature Ts°\ 

In order to investigate the dependence of the oscillations on the paramter A = L\/L^ we did a series 
of calculations using v = 3 and AM nuc = 0.001. In this case we find that oscillations occur only when 
the luminosity ratio A = Li/L^ lies in the range 0.84 < A < 0.92. We can also estimate the range of 
values of A that give oscillations from a linear stability analysis of the governing nonlinear equations (see the 
Appendix). In this case, the linear analysis predicts growing oscillations in the range 0.854 < A < 0.907. We 
have not found relaxation oscillations for A < 0.8, even though we have carried out a number of calculations 
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in this range. For example, with A = 0.3 and v = 10, the system goes through a single flash and then 
stabilizes at an HBRS luminosity lower than the initial assumed steady-state luminosity. Evidently this is 
not a self-consistent solution of the problem. Further, neglect of the luminosity L\ emerging from below 
the HBRS (i.e., setting A = 0) prevents matter in the envelope from re-heating when the luminosity drops 
and mass loss shuts down. Conversely, if the (presumed constant) interior luminosity L\ is too large, the 
system goes through a single thermonuclear pulse and then settles down into a state of steady burning. For 
example, with Li/L^ = 0.95, v = 3, Tf sm /r t h = 2.5, and r cx /r t h = 3000, there is no oscillation, and the shell 
temperature and luminosity stabilize at the assumed steady-state values. Thus, there appears to be a rather 
small region of parameter space in which the star exhibits such fuel-supply-limited relaxation oscillations. 

We have studied the dependence of the relaxation oscillations upon the assumed value of AM nuc , taking 
A = 0.9, v = 3, and s = 0.01 for all models. As shown in Table 3 below, we found that the oscillation 
period P decreases and the peak values of the temperature and luminosity decline as AM nuc decreases. The 
period appears to approach an asymptotic value ~ (15 — 17) x r t h for small HBRS masses. For the model in 
Table 3 with the lowest listed value of AM nuc , the amount of energy available from the shell flash is so small 
that the oscillations actually decay slowly in time, rather than attaining a roughly constant amplitude. 



Table 3. Variations with AM, 



AM nuc /M 


P/rth 






0.001 


33 


2.8 


1.36 


0.0006 


26 


2.1 


1.26 


0.0005 


22 


1.8 


1.20 


0.0004 


20 


1.6 


1.16 


0.0003333 


19 


1.5 


1.14 


0.0003077 


18 


1.25 


1.08 


0.0002857 


17-15 


: 


1.03 



We have also studied the dependence of the relaxation oscillations upon the assumed value of v, taking 
AM nuc = O.OO1M , A = 0.9, rj 0) = 10 7 K, and s = 0.01 for all models. As shown in Table 4 below, 
the system does not depend very much upon v, apparently because the temperature spikes so rapidly. We 
found no oscillations for v > 7. Within the limitations of our simplified model, a stronger dependence of 
the thermonuclear rate upon T apparently liberates energy so rapidly that the burning simply adjusts to a 
different stable-burning state. For calculations with v = 6, the peak temperatures reach ~ 5 x 10 8 K, and 
our approximation for the nuclear burning rate is invalid. 



Table 4. Variations with v 



V 


P/r th 




Lf ax /L { s 0) 


3 


33 


2.8 


1.36 


4 


35 


3.5 


1.44 


5 


38 


3.9 


1.47 


6 


51 


3.9 


1.48 



One aspect of the results surprised us: the conditions under which we found oscillations correspond 
to very low values of the nuclear-energy-generation exponent v. The values for which we find oscillations 
correspond to energy generation by the pp chain, whereas we had expected that the CNO cycle would 
dominate. Conceivably, this might also be a consequence of the terminal stage of shell burning. 

Finally, we note that in order to justify our suggested link between these relaxation oscillations and 
ejection of mass shells, the time scale for acceleration from rest to the escape velocity must be less than 
the oscillation period. This condition is easily satisfied; the acceleration time for material illuminated from 
below by luminosity £* ^ L E dd is ~ 0.13 year x (M*/M ) _1 / 2 (i?*/lO 13 cm) 3 / 2 (L»/i E dd - 1) _1 , assuming 
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the Thomson cross section as a conservative bound. This acceleration time is always much less than the 
^ 10 3 year oscillation period we find, so that we are well justified in applying our model to PN systems with 
evidence for periodic ejecta. 

6. Conclusions 

We have shown that fuel-supply-limited relaxation oscillations may develop in some stars near the ends 
of their nuclear-active lifetimes just prior to the ejection of a planetary nebula. In the benchmark model we 
discuss in this paper, with M* M©, L* ~ 10 4 _Lq, and \M\ ~ 1CP 5 Mq yr~ x , such oscillations occur when 
the luminosity generated by the H-burning relaxation shell (HBRS) falls to ~ 10% of the total luminosity 
of the star. The pulsation timescale we find to be ~ 72r t h ~ 57rf sm ~ 1400 years, closely comparable 
to the timescale inferred for the observed circumnebular shells (Balick et al. 2001; Zijlstra et al. 2002), 
for the present choice of parameters. While our choice of parameters does produce results consistent with 
the observed circumnebular shells, it is not clear whether such parameters are actually attained during the 
evolution of a real star. To determine whether such oscillations occur in more realistic stellar models, or 
whether they do indeed provide the correct explanation for the observed shells surrounding some planetary 
nebulae, will require much more detailed calculations of stellar evolution with mass loss than have been done 
to date. 
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Figure Caption 

Fig. 1. Fuel-supply-limited relaxation oscillations of the model described in §4 (model F in Tables 1 and 
2). The dotted curve represents the dimensionless temperature y = T s /Ts°\ the dashed curve gives the 
dimensionless HBRS luminosity L s /L^\ and the solid curve is the dimensionless H abundance x = X s /Xs°^ 
in the HBRS. The abscissa gives the time elapsed from the start of the calculation in units of the thermal 
timescale defined by equation (24), t/rth- At the onset of a pulse, the rapid increase in luminosity drives 
the H abundance in the HBRS to zero in less than ~ 39rth ~ 780 years. During this initial phase, the shell 
temperature rises rapidly to a peak value ~ 62 times the intial value. The energy generated during this 
pulse is gradually radiated away during the next ~ 39T t h ~ 780 years, and the luminosity declines toward 
the initial luminosity on the same timescale. When the surface luminosity and the associated mass-loss rate 
have declined sufficiently, the H abundance in the HBRS begins to build up again. In ~ 72r t h ~ 1400 years 
<~ 57rf sm , the energy generation rate in the HBRS builds up to a level exceeding the level of steady-state 
burning, and another pulse occurs. 

Appendix: Linear Stability Analysis 

To investigate the stability of the equilibrium solution x = y = 1 of equations (38), let x = 1 + £ and 
y = 1 + rj and expand the functions g(y) and h(y) in Taylor series around y = 1 in the form 

.9(1 + 17) = l + b v + 0( V 2 ), h(l + r}) = l + c^ + O^ 2 ), (Al) 
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where the coefficient b depends on the parameter sq and the coefficient c depends on the parameters sq and 
v. Substituting these expressions into equations (38) and linearizing yields the equations 

^- = - (7i & + 72c)»7, {Ala) 

f| = (l_A)£+[c(l-A)-%. (A26) 

If we assume a nonzero solution of these equations in the form £ = £oexp(ai), r\ — rjoexp(oct), then we find 
that a satisfies a quadratic equation with roots 

a = [B ± (B 2 - 4C) 1/2 ]/2, (A3) 

where 



B = (1- A)c-6-72, C=[(l-A)7i+7 2 ]6. (A4) 

The solution is an oscillation when £> 2 — 4C < and these oscillations are growing when B > 0. 

Using this analysis, we can determine the regions of parameter space in which we expect to find growing 
oscillations. For example, with so = 0.03, v = 3, 71 = 0.8, and 72 = 0.001143 (as in model F in Table 1), we 
find that the solutions of the linearized equations are oscillatory for values of the parameter A in the range 
0.807 < A < 0.957, and these oscillations are growing for A in the range 0.807 < A < 0.908. 
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